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Abstract 

Many uncertainty sets encountered in control systems analysis and design can be 
expressed in terms of semialgebraic sets, that is as the intersection of sets described 
by means of polynomial inequalities. Important examples are for instance the solu¬ 
tion set of linear matrix inequalities or the Schur/Hurwitz stability domains. These 
sets often have very complicated shapes (non-convex, and even non-connected), 
which renders very difficult their manipulation. It is therefore of considerable im¬ 
portance to find simple-enough approximations of these sets, able to capture their 
main characteristics while maintaining a low level of complexity. For these reasons, 
in the past years several convex approximations, based for instance on hyperrect¬ 
angles, polytopes, or ellipsoids have been proposed. 

In this work, we move a step further, and propose possibly non-convex approx¬ 
imations, based on a small volume polynomial superlevel set of a single positive 
polynomial of given degree. We show how these sets can be easily approximated 
by minimizing the norm of the polynomial over the semialgebraic set, subject 
to positivity constraints. Intuitively, this corresponds to the trace minimization 
heuristic commonly encounter in minimum volume ellipsoid problems. From a com¬ 
putational viewpoint, we design a hierarchy of linear matrix inequality problems to 
generate these approximations, and we provide theoretically rigorous convergence 
results, in the sense that the hierarchy of outer approximations converges in volume 
(or, equivalently, almost everywhere and almost uniformly) to the original set. 

Two main applications of the proposed approach are considered. The first one 
aims at reconstruct ion/approximation of sets from a finite number of samples. In 
the second one, we show how the concept of polynomial superlevel set can be used to 
generate samples uniformly distributed on a given semialgebraic set. The efficiency 
of the proposed approach is demonstrated by different numerical examples. 
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1 Introduction 


In this paper, we address the problem of how to determine “simple” approximations of 
semialgebraic sets in Euclidean space, and we show how these approximations can be 
exploited to address several problems of interest in systems and control. To be more 
precise, given a set 

/C = {a: e M"" : gi{x) >0, i = 1, 2,..., m} (1) 

which is compact, with non-empty interior and described by given real multivariate poly¬ 
nomials gi{x),i = 1,2,... ,m, and a compact set B D 1C, we aim at determining a so-called 
polynomial superlevel set (PSS) 

U{p) = {x ^ B -. p{x) > 1}. (2) 

that constitutes a good outer approximation of the set 1C of interest and converges strongly 
to K, when increasing the degree of the real multivariate polynomial p to be found. 

In particular, the proposed PSS is based on an easily computable polynomial approxi¬ 
mation of the indicator function of the set 1C. In the paper, we show that suitable ap¬ 
proximations of the indicator function can be obtained by solving a convex optimization 
problem whose constraints are linear matrix inequalities (LMIs) and that, as the degree 
of the approximation increases, one converges in L^-norm, almost uniformly and almost 
everywhere to the indicator function of the semialgebraic set 1C of interest. Moreover, the 
set approximations provided in this paper can be thought as a direct generalization of 
classical ellipsoidal set approximations, in the sense that if second degree approximations 
are used, we exactly recover well-known approaches. 

The main motivation for the problem addressed in the paper is the fact that semialge¬ 
braic sets are frequently encountered in control. As an example, consider the Hurwitz or 
Schur stability regions of a polynomial. It is a well-known fact that the these regions are 
semialgebraic sets in the coefficient space. The polynomial inequalities that define these 
stability sets can be derived from well-known algebraic stability criteria. Another classical 
example of semialgebraic sets arising in control are LMI feasibility sets, also called spec- 
trahedra. Indeed, LMI sets are (convex) basic semialgebraic sets. To see this, consider 
the LMI set 

^LMi = { 3 ; G M” : F{x) = Fo + FiXi -|- • • • -|- F^Xn F 0} 

where the matrix F{x) has size m x m, and observe that a vector x belongs to /Clmi if 
and only if all the coefficients of the univariate polynomial 

s I—>■ det {sIm + F(x)) = gi{x) g 2 {x)s -|- • • • -I- gjn{x)s^ ^ s^ 

are nonnegative, i.e. x belongs to the set K is defined in 0. where the polynomials gi{x) 
are by construction sums of principal minors of the matrix F[x). The approach taken 
in this paper is the following: given the set /C, we search for a minimum volume PSS 
that contains the set /C. Since there is in general no analytic formula for the volume of 
a semialgebraic set, in terms of the coefficients of the polynomials defining the selQ it is 

^See however reference [33] which explains how explicit formulas can be obtained with discriminants 
in exceptional cases. 
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very challenging to solve this optimization problem locally, let alone globally. Instead, the 
main contribution of this paper is to describe and justify analytically and geometrically 
a computationally tractable heuristic based on L^-norm or trace minimization. Second, 
we show that the same approach can be employed to obtain the largest (in terms of the 
L} surrogate for the volume) PSS inscribed in /C. Moreover, it is shown how the ideas 
put forth in this paper can be used to address two important problems: i) reconstruc¬ 
tion/approximation of a (possibly non-semialgebraic) set from samples belonging to it, 
and ii) uniform generation of samples distributed over a semialgebraic set. Examples of 
applications in a systems analysis and controller design context are also provided. 

The work presented in this paper is an extension of the preliminary results in the con¬ 
ference papers m and [12] , and it provides a more in depth analysis of both theoretical 
and implementation aspects. In particular, with respect to [H], the present manuscript 
contains more detailed proofs of the theoretical results, provides detailed algorithmic de¬ 
scriptions, and introduces inner PSS approximations. Similarly, the results on random 
sample generation of na are here described in more details, and an algorithm is pro¬ 
vided. Finally, all examples in the paper are new, and more control oriented applications 
are considered. 


1.1 Previous work and related literature 

The idea of approximating overly complicated sets by introducing simpler and easy man¬ 
ageable geometrical shapes is surely not new, it has a very long history, and it arises in 
different research fields such as optimization, system identification and control. In partic¬ 
ular, in the systems and control community, the most common approach is to introduce 
outer bounding sets, that is sets of minimum size which are guaranteed to contain the set 
to be approximated. For instance, in the context of robust filtering, set-theoretic state 
estimators for uncertain nonlinear dynamic systems have been proposed in [Il[l8l[20l|39|. 
These strategies adopt a set-membership approach [T9l |38| , and construct (the smallest) 
compact set guaranteed to bound the system states that are consistent with the measured 
output and the norm-bounded uncertainty. The most common geometrical shape adopted 
in these work is the ellipsoidal one, for the double reason that it has a very simple descrip¬ 
tion - the center and the shape matrix are sufficient to provide a complete characterization 
- and that its determination usually can be formulated as a convex (usually quadratic) 
optimization problem. The use of ellipsoidal sets in the state estimation problems was 
introduced in the pioneering work [HH| and used by many different authors from then on; 
see, for example, [lEl ED] • Outer approximation also arise in the context of robust fault 
detection problems (e.g., see [26|) and of reachability analysis of nonlinear and/or hybrid 
systems [251128] . Similarly, inner approximations are employed in nonlinear progra mming 
[M] , in the solution of design centering problems m and for fixed-order controller design 
[21]. In this case, one aims at constructing the set largest size inscribed in the set of 
interest. 

Besides ellipsoids, other shapes have been considered in the recent literature. The use 
of polyhedrons was proposed in [27] to obtain an increased estimation accuracy, while 
zonotopes have been also recently studied in PEI]. In P a heuristic based on polynomial 
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optimization and convex relaxations is proposed for computing small volume polytopic 
outer approximations of a compact semialgebraic set. More recent works, like for instance 
[SlEaE2], employ sets defined by semialgebraic conditions. The closest approach to the 
one proposed in our paper can be found in |H2], in which the authors use polynomial 
sum-of-squares (SOS) programming to address the problem of fitting given data with a 
convex polynomial, seen as a natural extension of quadratic polynomials and ellipsoids. 
Convexity of the polynomial is ensured by enforcing that its Hessian is matrix SOS, and 
volume minimization is indirectly enforced by increasing the curvature of the polynomial. 
In [S] the authors propose moment-based relaxations for the separation and covering 
problems with semialgebraic sets, thereby also extending the classical ellipsoidal sets used 
in data fitting problems. 

Recently, the authors of [Hj have proposed an approach based on randomization, which 
constructs convex approximations of generic nonconvex sets which are neither inner nor 
outer, but they enjoy some specific probabilistic properties. In this context, an approx¬ 
imation is considered to be reliable if it contains “most” of the points in the given set 
with prescribed high probability. The key tool in this framework is the generation of ran¬ 
dom samples inside the given set, and the construction of a convex set containing these 
samples. 


1.2 The sequel 

The paper is organized as follows. In Section the notation used in this paper is intro¬ 
duced and the central problem addressed in this paper is defined. In order to be able 
to numerically solve the set approximation problem of interest, in Section a related 
polynomial optimization problem is introduced and numerical methods for solving it are 
described in Section In Section we discuss how the results in this paper can be 
used to find inner approximations of semialgebraic sets. A first set of numerical exam¬ 
ples is provided in Section]^ Using the central results on set approximation mentioned 
above, in Section we address the problem of reconstructing a set from a finite number 
of points in its interior. In Section we provide algorithms for uniform sample generation 
in semialgebraic sets and in Section some closing remarks are provided. 


2 Problem statement 

Before a description of the main problem addressed is provided, we introduce the basic 
notation that is used throughout the paper. 

2.1 Notation 

The notation A >- 0 0) means that the symmetric matrix A is positive definite (semidef- 

inite), and given two matrices A and B we write A ^ B whenever A — B ^ 0. Given a 
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set /C C M”, its indicator function is defined as 




1 if a: e /C 

0 if X ^ K, 


(3) 


and its volume or, more precisely, the Lebesgue measure of /C, is denoted by 


vol 1C = 



Iic{x)dx. 


The set of all real coefficient polynomials of degree less than or equal to d is denoted by 
Pd- The monomial basis for this set is represented by the (column) vector e Pd, so 
that any p & Pd can be expressed in the following form 


P{^) = ix)p = 7rf^/2] (x)P7r rd /21 (a:) 

where p is a real (column) vectoi|^ and P is a symmetric matrix of appropriate size, often 
referred to as Gram matrix. Also, we denote by Ti 2 d the set of polynomials p G P 2 d that 
can be represented as sums of squares of other polynomials, i.e. 


Up 

P = '^pt Pk^ Pd, k = l,...,np. 

k=l 

Finally, given a polynomial p G Pd, define its norm over a compact set B, denoted by 
Z/g or just when the set B used is clear from the context, as 


Iblli 


p{x)dx. 


2.2 Problem Statement 

With the notation defined above, we are now ready to define the central problem in this 
paper. We consider the basic semialgebraic set JC defined in Q, which is assumed to be 
compact and with a non-empty interior. 

As discussed in the Introduction, the set K, has typically a complex description in terms 
of its defining polynomials (e.g. coming from physical measurements and/or estimations). 
For this reason, we aim at finding a “simpler” approximation of this set which has enough 
degrees of freedom to capture its characteristics. This approximation is the polynomial 
superlevel set (PSS) U{p) defined in ([^ in terms of a real multivariate polynomial p E Pd 
of given degree d. This degree controls the complexity of the approximation. Among the 
family of possible PSS that can be constructed, we search for the one that provides the 
set V({p) of minimum volume while containing the set of interest /C, hence capturing most 
of its the geometric features. Formally, we define the following optimization problem 

^Note that we use p to denote both the polynomial and the vector of its coefficients whenever no 
ambiguity is possible. 
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Problem 1 (Minimum volume outer PSS) Given d & N and a compact semialge- 
braic setJC, find a polynomial p G Pd whose PSSU{p) is of minimum volume and contains 
JC. That is, solve the following optimization problem 

v*^ = inf Yo\U{p) 

pePd (4) 

s.t. /C C U{p). 

Note that this problem can be viewed as the natural extension of the problem of computing 
the minimum volume ellipsoid containing 1C. Indeed, if JC is convex and the polynomial 
p is quadratic (d = 2), then the infimum of problem Q is attained, and the optimal set 
U{p) is given by the unique (convex) ellipsoid of minimum volume that contains /C, called 
Lowner-John ellipsoid. In particular, if K, is the convex-hull of a finite set of points, this 
ellipsoid can be computed by convex optimization, see e.g. [H §4.9]. 

We remark however that, for d greater than 2, the optimization problem Q is nonlinear 
and semi-infinite, in the sense that the optimization is over the finite-dimensional vector 
space Pd, but subject to an infinite number of constraints, necessary to cope with the set 
inclusion. 

Theorem 1 The sequence of infima of problem 0 monotically converges from above to 
vol K, i.e. for all d> 1 it holds and hm(i^.oo = vol 1C. 

Proof: As in |23l Section 3.2], let x ^ d(x,)C) be the Euclidean distance to set 1C and 
with Cfc > 0 let ICe := {x & B : d{x,lC) < e^} be an open bounded outer approximation 
of 1C, so that B\lCk is closed with hmfc^.oo = 0. By Urysohn’s Lemma [371 Section 12.1] 
there is a sequence of continuous functions {fk)kef^ with fk '. B ^ [0,1] such that fk = 0 
on B\lCk and /^ = 1 on 1C. In particular, notice that vol/C < vo\U{fk) < vollC + vo\lCk\lC 
and since limfc_j.oo vol lCk\lC = 0 it holds limfc_j.oo volZ//(/fc) = vol 1C. 

By the Stone-Weierstrass Theorem m Section 12.3] we can approximate fk uniformly on 
B by a sequence of polynomials {p'kJdGN with G Pd, i.e. sup^g^ \fk{x) - p'k^fx)\ < e'^ 
with limd^oo4 = 0. Defining pk^ ■= p'k,d + 2ed, the sequence of polynomials {pk,d)deN 
converges uniformly to fk from above, i.e. pk^ > fk on B and limc;_j.oo sup^gg |/fc(a:) — 
Pk,d{^)\ = 0- This implies that := volU{pk,d) > yolU{fk) and limd^oovl^ = yo\U{fk). 
Recalling lim^^oo vol W(/fc) = vol/C, it follows that linik^^ooPk d ~ vol/C which proves, 
up to extracting a subsequence indexed by d, the existence of a minimizing sequence of 
polynomials for optimization problem Q. 

Finally, the inequality readily follows from the inclusion Pd C Pd+i- D- 

Before introducing the approach we propose for the solution of Problem [T| in the next 
subjection we briefly recall some recent results which are closely related to the problem 
considered in this paper, for the special case of homogeneous polynomials. 


2.3 Remark on a convex conic formulation 

In this section, we summarize existing results for the case when the polynomial q £ Pd 
defined as q{x) = 2 — p{x) is assumed to be a homogeneous polynomial, or form, of even 
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degree d = 2S in n variables. First note that, with this change of notatiorj^ the PSS U{p) 
corresponds to the unit sublevel set V{q) ={x G M” : q{x) < 1} of the polynomial q. In 
[HOI Lemma 2.4] it is proved that, when q is homogeneous, the volume function 

q I—7- vol V{q) 

is convex in q. The proof of this statement relies on the striking observation [H3] that 

vol V{q) = Cd [ e-^^^^dx 

jRn 

where Cd is a constant depending only on d. Note also that boundedness of V{q) implies 
that q is nonnegative, since if there is a point xq G such that q{xQ) < 0, and hence 
xq G V{q), then by homogeneity of q it follows that q{Xxo) = \‘^^q{xQ) < 0 for all A and 
hence Axq G V{p) for all A which contradicts boundedness of V{p). This implies that 
problem Q, once restricted to nonnegative forms, is a convex optimization problem. 

Moreover, in m Lemma 2.4] explicit expressions are given for the first and second order 
derivatives of the volume function, in terms of the moments 

[ x^e-^^^^dx (5) 

for a G |qi| < 2d. In an iterative algorithm solving convex problem 0 , one should 
then be able to compute repeatedly and quickly integrals of this kind, arguably a difficult 
task. Moreover, when q is not homogeneous, we do not know under which conditions on 
q the function vol V{q) is convex in q. 

Motivated by these considerations, in the remainder of this paper we propose a simpler 
approach to the solution problem 0 , which is not restricted to forms, and which does 
not require the potentially intricate numerical computation of moments (j^ of exponen¬ 
tials of homogeneous polynomials. The introduction of this approach is motivated by its 
analogy with the well-known trace heuristic for ellipsoidal approximation, and it consists 
of approximating the volume by means of the L^-norm of the polynomial p. 


3 L^-norm minimization 

It is assumed that a “simple set” B C containing /C is known. By “simple” we 
mean that analytic expressions of the moments of the Lebesgue measure on B should be 
available, so that integration of polynomials can be carried out readily. In the following, 
we assume that the set B is an n-dimensional hyperrectangle of the form 

B = [a, b] = {x E < Xi < bi, i = 1,2,, n} (6) 

with a and b given vectors of E". This is a very mild assumption since, given a semial- 
gebraic set like the set K, above, one can easily compute an hyperrectangle containing it; 

^ The polynomial g = 2 — p is introduced because the results in m are derived for sublevel sets, not 
super level sets. 
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see Section |4.1| for details. We note that more complex sets B ^ K, can be considered, 
provided that integration of polynomials over it is easily done. 


Assnme now, without loss of generality, that the polynomial p used to build the PSS 
is non-negative on B. Then, observe that by definition of PSS (see Figure for an 
illustration) we have 

P > h{p) on B. 



Figure 1: Illustration of Chebychev’s inequality: the polynomial is always greater or equal 
than the indicator function of p{x) > 1, hence the integral of p over B is always an upper 
bound of the volume of U {p). 


Hence, integrating both sides we get the following inequality 


p{x)dx> / lu{p){x)dx = Yo\U{p). 

? Jb 


(7) 


This inequality is indeed widely used in probability, where it goes under the name of 
Chebyshev’s inequality, see e.g. m §2.4.9]. Note that, since the polynomial p is nonneg¬ 
ative on B, then the left-hand side of inequality Q corresponds to the L^-norm of p on 
B, so that the inequality simply becomes 


Iblli > voiz^(p). 


( 8 ) 


These derivations motivate us to the formulation of the following L^-norm minimization 
problem, which we choose as a surrogate of the original minimum volume outer PSS 
introduced in Problem 1. 


Problem 2 (Minimum L^-norm outer PSS) Given a semialgehraic set K,, a bound¬ 
ing set B ^ K, and a degree d, solve the optimization problem 

K = 

p&Pd 

s.t. p > 0 on B (9) 

p > 1 on /C. 


Note that a L^-norm minimization approach was originally proposed in [22] for the nu¬ 
merical computation of the volume and of the higher order moments of a semialgehraic 
set. The intuition underlying the formulation of Problem is similar. We now elaborate 
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on some of the characteristics of the minimum L^-norm outer PSS problem dehned above. 
First note that, for fixed d, when solving Problem we are minimizing an upper-bound 
on the volume of the PSS. Thus, the solution is expected to be a good approximation 
of the set /C. Second, it can be shown that, as the degree d increases, the Chebyshev 
bound it becomes increasingly tight. Indeed, the following fundamental result shows 
that the proposed solution converges to the minimum volume outer PSS. 

Theorem 2 Given d G N, the infimum in problem I® is attained for a polynomial p*^ G 
Pd- Moreover, w*^ > v*^ andUijff) 3 /C. Finally andlim.d^oo'Wd = ~ 

vol /C. 


Proof: Let us first extend optimization problem Q to continuous functions: 


w* = inf / f{x)dx 

f Jb 

s.t. feC+{B) 

/-lGC+(/C) 


( 10 ) 


where C+{B) denotes the convex cone of non-negative continuous functions on B. Observe 
that since / is non-negative on B, the objective function ||/||i = J^f{x)dx is linear. 


= sup 

s.t. 


Problem (10) is an infinite-dimensional linear programming (LP) problem in cones of 
non-negative continuous functions. It has a dual LP, in infinite-dimensional dual cones of 
measures: 

J ii{dx) 

fi(dx) + jl(dx) = Is{x)dx (11) 

A e C'^B) 

h e c;(/C) 

where Cf{B) is the cone of non-negative continuous linear functionals on C+{B), identified 
with the cone of Borel regular non-negative measures on B, according to a Riesz Repre¬ 
sentation Theorem [HU Section 21.5]. In LP ( [IT| ) the right hand side in the equation is the 
Lebesgue measure on B. Since the mass of non-negative measures fi and fi is bounded, 
it follows from Alaoglu’s Theorem on weak-star compactness m Section 15.1] that the 
supremum is attained in dual LP ( [IT| ) and that there is no duality gap between the primal 
and dual LP, i.e. v* = w*, see also e.g. [SI Theorem IV.7. 2]. 

Moreover, as in the proof of [2S1 Theorem 3.1], it holds v* = vol /C. To see this, notice 
first that the constraint fi + fi = Ib jointly with pL G F(,_(/C) imply that ja < Ik. and hence 
J /i < J Ik = vol JC for every /i feasible in LP In particular, this is true for an 

optimal Id* attaining the supremum, showing J /a* = v* < vol /C. Conversely, the choice 
fi = Ik is trivially feasible for LP (11) and hence suboptimal, showing v* > f fa = vol/C. 
From this proof it also follows that the only optimal solution to LP ( [TT| is the pair 
{ft*, ft*) = (I/Cilex/c)- 

Now let us prove the statements of the Theorem: 

• Attainment of the infimum in problem ([^ follows from continuity (actually linearity) 
of the objective function ||p||i = f^p(x)dx = 0 which is a norm (i.e. ||p||i = 0 implies 
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p = 0 for p G Pd) and compactness of the set {p G : p G C{B), ||p||i < r} for any 
fixed r > 0. 


• uj'd — '^d follows readily from Q. 

• w*^> follows readily from Pd C Pd+i- 


• Finally, limd-j.oo = vol/C is a consequence of v* = w* = vol/C (proven above) and 
the Stone-Weierstrass Theorem [371 Section 12.3] allowing to approximate uniformly 
on B by polynomials any continuous function in a minimizing sequence for LP (10), 
i. e. lim^^oo = w*. 


□. 

Some remarks are at hand regarding the above result, which represents one of the main 
contributions of the paper. 


Remark 1 (Convergence almost everywhere) Note that Theorem^implies that, for 
high enough order of approximation, the PSS obtained by minimizing the L^-norm of the 
polynomial defining it ean be “arbitrarily elose” to the semialgebraie set of interest. More 
preeisely, as d ^ oo, ||p^||i and, as a consequenee volW(p^), eonverges to vol/C. Sinee, 
K, C U{p*fi), the Lebesgue measure of the difference between these sets converges to zero. 
In other words, one has almost everywhere convergence. From Theorems 2.5.1 and 2.5.3 
in m the convergence is also almost uniform, up to extracting a subsequence. 


Remark 2 (Trace minimization) We provide a geometric interpretation that further 
justifies the approximation of the minimum-volume PSS with the minimum L^-norm PSS. 
To this end, we first note that the objective function in (§ reads 



p{x)dx 


7rJ(x)PTTs{x)dx = trace ( P / 7rs{x)7rJ(x)dx 


trace PM (12) 


where 


M 


= / tts{x)tt^ 

Jb 


{x)dx 


is the matrix of moments of the Lebesgue measure on B in the basis 7r^(a:). Note that, if 
the basis in equation (12) is chosen such that its entries are orthonormal with respect to 
the (scalar product induced by the) Lebesgue measure on B, then M is the identity matrix 
and inequality becomes 

trace P > yo\U{p) 


which indicates that, under the above constraints, minimizing the trace of the Gram matrix 
P entails minimizing the volume ofU{p). It is important to remark that, in the case of 
quadratic polynomials, i.e. d = 2, we retrieve the elassical trace heuristic used for volume 
minimization of ellipsoids, see e.g. m- Indeed, if B = [—1,1]"', then the basis 'rfix) = 
^x is orthonormal with respeet to the Lebesgue measure on B and ||p||i = |trace P. 
Moreover, note that the eonstraint that p is nonnegative on B implies that the eurvature 
of the boundary ofU{p) is nonnegative, henee thatU{p) is eonvex. Thus, U{p) is indeed 
an ellipsoid. 
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Remark 3 (Choice of B) We finally remark that, as previously noted, the assumption 
of the bounding set B being an hyperrectangle can be easily relaxed. Indeed, in order to 
develop a computationally manageable optimization in Problem 2, B can be selected as a 
semialgebraic set, provided that the polynomials defining the set should be such that the 
objective function in problem © is easy to compute. In particular, if 

P{x) = ^dix)P = '^PaMx)]a 


then 

p{x)dx = '^Pa [Trdix)]adx = ^PaVa 

a a 

and we should be able to compute easily the moments J^[TTdix)]adx of the Lebesgue measure 
on B with respect to the basis Udix). 



4 LMI hierarchy to compute the PSS 


In this section, we provide the basic details on the numerical computation of the solution 
of the minimum L^-norm PSS introduced in Problem 2. Note that, in problem ([^, we aim 
at finding a polynomial p & Pd such that i) p is positive on B, and ii) p — 1 is positive on 
/C. In order to obtain a numerically solvable problem, we enforce positivity by requiring 
the polynomial to be SOS, and use Putinar’s Positivstellensatz; e.g., see [Ml ESI ISl 1^ - 
More precisely, fix r G N, and consider the problem 


W^r d = Him 
pePd. 


p{x)dx 


s.t. 


p{x) = So,b{x) + Sj^B{x){xj - aj){bj 
So,B £ ^2r 

Sj,B £ S2(r-1), i = 2, . . . , n 

m 

p{x) - 1 = so,x;(a:) + Si^!c{x)gi{x) 

i=l 

So,IC G S2r 

^i,IC £ ^2(i—rp) I 1; 2, , 771. 


(13) 



p{x) positive on B= [a, b] 


> p{x) — 1 positive on K, 


where r* is the smallest integer greater than half the degree of pj for i = 1, 2,..., m. It 
should be noted that the objective function of problem (13) is an easily computable linear 
function of the coefficients of the polynomial p. Moreover, the constraints can be recast 
in terms of Linear Matrix Inequalities (LMIs); see, for instance, (SHI)- Several Matlab 
toolboxes have efficient and easy to use interfaces to model problems of the form above; 
e.g., see YALMIP [3T]. 
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Not only we can numerically solve problem (13), but the following result holds. This 
theorem is an immediate consequence of the results in [HHj . 


Theorem 3 Let us denote hyp^rd ® solution of problem (13). Then, the following hold 


i) for each d G N, the value of problem (13) converges to the value of problem 0 as 
r ^ oo, i.e. lim^^oo''^ 2 r,d = '^d> 

a) for any 2r > d, p^^d ^0 on B, 

Hi) for any 2r > d, P 2 rd on K. 

We conclude that be used to compute a PSS approximation for /C. For our 

numerical examples, we have used the YALMIP EH interface for Matlab to model the 


LMI optimization problem (13) and the SDP solver SeDuMi mi to numerically solve the 
problem. Since the degrees of the semialgebraic sets we compute are typically low (say 
less than 20), we did not attempt to use alternative polynomial bases (e.g. Chebyshev 
polynomials) to improve the quality and resolution of the optimization problems; see [23] 
for a discussion on these numerical matters in the context of semialgebraic set volume 
approximation. 


4.1 Computing Bounding Box B 

As noted in jS] Remark 1], an outer-bounding hyper-rectangle B = [a,b] of a given semial¬ 
gebraic set K can be found by solving relaxations of the following polynomial optimization 
problems 

Oj = arg min Xj subject to a: G /C, j = 1,..., n, 
bn = argmaxxo subject to a: G /C, j = 1, ...,n, 

a;GR’» 

which compute the minimum and maximum value of each component of the vector x over 
the semialgebraic set /C. 

To illustrate how this can be done, let us concentrate on approximating the value of 
Oj. First, note that the problem of computing aj is equivalent to solving the following 
polynomial optimization problem 

Oj = maxy subject io Xj — y >t) for all x & 1C. 

Then, formulate the following convex optimization problem 

aj^2r = max y 

m 

s.t. Xj-y = So(x) -F Si{x)gi{x) 

1=1 

50 £ Tj2r’i 

51 G Tj2(r—ri)i i 1)2,..., 171. 

Using the same reasoning as above, it can be shown that: i) aj^ 2 r < CLj for ah r, and ii) 
limr^.oo Oj, 2 r = (ij- Moreover, the problem above can be recast as an LMI optimization 
problem. 
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5 Inner approximations 


The approach described in the previous sections can be readily extended to derive inner 
approximations of the set /C, in the spirit of [ElEaEl!. The idea is just to construct an 
optimal outer PSS of the complement set 

K. = B\)C = {x eR"' : gi(x) < 0 or • • • or gmix) <0, i = 1,2,..., m} fi B 


= (/Cl u ;C 2 u • • • u /c^) n b. 


with )Cj = {a: e M” : gj{x) < 0}. 

Note that, since the set whose indicator function we want to approximate is a union of 
basic semialgebraic sets, the optimization problem to be solved becomes 


min 

iibiii 



p^Pd 






s.t. 

p 

> 

0 

on 

B 


p 

> 

1 

on 

/Cl 


p 

> 

1 

on 

/C 2 


p 

> 

1 

on 

/C™ 


(14) 


and let attain the minimum. The corresponding optimal inner approximation is given 
by the polynomial sublevel set 


V(Pd) = {xe B:pa{x) < 1}. 


In this case, one can think of the polynomial 1 — as a lower bound for the indicator 
function of the set 1C. 


Given the fact that the optimization problem (14) provides an outer approximation of the 
set /C, one has the following result whose proof is similar to that of Theorem]^ 


Corollary 1 For all d eN it holds V(p^) C /C. Moreover lim^^oo vol V (p^) = vol /C. 


As before, to be able to numerically approximate the solution of problem 
place its polynomial positivity constraints by their LMI approximations as 
Section HI 


(14), we re¬ 
described in 


6 Numerical examples 

In this section, we present several examples that illustrate the performance of the proposed 
approach. 
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6.1 Discrete-time stabilizability region 


As a control-oriented illustration of the PSS approximation described in this paper, 
consider m Example 4.4] which is a degree 4 discrete-time polynomial z e C i—>■ 
X 2 + 2xiZ — {2xi + X 2 )z^ + to be stabilized by means of 2 real control parameters 
Xi,X 2 . In other words, we are interested in approximating the set K, of values of X\,X 2 
such that this polynomial has its roots with modulus less than one. An explicit basic 
semialgebraic description of the stabilizability region is built using the Schur stability 
criterion, resulting in the following basic semialgebraic set: 


1C = {a: e : 5fi(x) = 1-|-2x2 > 0, (15) 

92 { x ) = 2 — 4xi — 3x2 > 0 , 

5-3(x) = 10 — 28 xi — 5x2 — 24 xiX 2 — 18x2 > 0, 

g ^ ix ) = 1 — X2 — 8x^ — 2 xiX 2 — X2 — 8x^x2 — 6x1X2 > 0}. 

This set is nonconvex and it is included in the box B = [—0.8, 0.6] x [—0.5,1.0]. In Figure 
we represent the PSS outer approximations of /C for d = 6 and d = 12 respectively, 
while Figure shows the graph of the degree d = 12 polynomial Pi 2 ,i 2 (^) constructed by 
solving optimization problem (13) with 2r = d. 


As discussed before, we can also use the approach proposed in this paper to obtain inner 
approximations of /C. In Figure we depict the inner approximation obtained using 
optimization problem (14) with 2r = d = 8. 




Figure 2: Degree 6 and degree 20 outer PSS approximation (red) of stabilizability region 
1C (inner surface in light blue). The green box corresponds to the bounding set B. 
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Figure 3: Degree 20 polynomial approximation (upper surface in red) of the indicator 
function (lower surface in blue) of the nonconvex planar stabilizability region /C. 

6.2 PID stabilizability region 

We now turn our attention to an example related to fixed order controller design. Consider 
0 Example 2.2], in which the authors examine the problem of stabilizing the plant 
Pis) = ^ where 

Nis) = — s — 1; 

Dis) = s® + + 32s^ + 26s^ + 65s^ - 8s + 1. 

by means of a PID controller of the form if’piD(s) = fcp + ^ + kos. In particular, they are 
interested in finding the set of stabilizing PID gains, that is the set of gains for which the 
closed-loop characteristic polynomial sD(s) -|- {ki -|- kps + kDs‘^)Nis) is Hurwitz. For this 
special class of controllers, the authors provide a method based on the so-called signature 
of a set of properly constructed polynomials to determine the set of all PID gains that 
stabilize the plant. One should note that this procedure is not easily generalizable to 
more general classes of fixed order controllers. 

In our setup, we are interested in approximating the set 

/C = {x G : sDis)+{ki+kps+kDs‘^)N{s) is Hurwitz, ki = 25(xi —1), kp = 10(x2—1.5), kp, = 10(x3- 
with bounding box B = [—1,1]^. As one can see in Figurethe approached proposed in 
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Figure 4: Left: degree 8 inner PSS approximation (red) of stabilizability region JC (inner 
snrface in light blue). Right: degree 8 polynomial approximation (upper surface in red) 
of the indicator function (lower surface in blue of /C. 

this paper provides a very good approximation of the set of stabilizing gains, even for a 
PSS of relatively low order {d = 14). 


7 Reconstructing/approximating Sets from a finite 
number of samples 

A particularly interesting case is when the semialgebraic set K, is discrete, that is, it 
consists of the union of N points 

N 

)C = C M". 

i=l 

This sitnation arises for instance when the objective is to try to approximate a given set 
(possibly non-semialgebraic) from a given number of points in its interior. An example 
of this is the reconstruction of reachable sets by using randomly generated trajectories. 
This setup is discussed in dlllSlEH]. 

From a computation viewpoint, an important feature is that, in the case of a discrete set, 
the inclusion constraint K QU{p) is equivalent to a hnite number of inequalities 

>1, i = 1,... ,N 

which are linear in the coefficients of p. This fact allows to deal with problems with 
rather large N. Moreover, in this latter case, where the number of points N is large while 
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Figure 5: Left: set of stabilizing PID gains. Right: its degree 14 optimal outer PSS 
approximation. 



the dimension n is relatively small, the constraint that p is nonnegative on B can also be 
(approximately) handled by linear inequalities 

enforced at a dense grid of points G B, for M snfhciently large. Hence, in this case one 
can construct a pure linear programming (LP) approach. Note that, even if this approach 
does not guarantee that p is nonnegative everywhere on it still ensures that K QU{p), 
which is what matters primarily in our approach. 

To illustrate the performance of the proposed method, we first consider N = 100 points 
in the box B = [—1, 1]^. The points are generated mapping Gaussian points with vari¬ 
ance 0.1/ and mean value chosen with eqnal probability between [0.4, 0.3]^, [—03, —0.5]^, 
[—0.5, 0.4]^. On Fignre|^we represent the solntions p of degrees 2, 5, and 9 of minimiza¬ 
tion problem Q. A few comments about the obtained solution are at hand. First, we see 
that the solution for d = 2 corresponds to the Lowner-John ellipsoid, see e.g. n §4.9]. 
Second, it can be observed that, as the degree of p increases, the set Vl{p) becomes dis¬ 
connected, so as to better capture the different regions where the points are concentrated. 
We note that, in the case of discrete points, it is not advisable to select high valnes of d, 
since indeed, in the limit, the optimal polynomial would correspond to a function with 
spikes corresponding to the location of the considered points. Finally, we remark that the 
possible side effects near the border of B on the right hand side figure can be removed by 
enlarging the bonnding set B. 

As a second illustrative example, we consider = 10 points inB = [—1,1]^. The solutions 
p of degrees 4, 6, 9, and 14 of minimization problem Q is depicted in Figure Here 
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Figure 6: Minimum L^-norm PSS at 100 points (blue), for degree 2 (left), 5 (center), and 
9 (right). 


too we observe that increasing the degree of p allows to capture point clusters in distinct 
connected components. 


8 Uniform sampling over semialgebraic sets 

In this section, we consider a problem that can be seen as the “dual” of the one considered 
in the previous section; that is, instead of trying to reconstruct/approximate the indicator 
function of an unknown set from points belonging to its interior, we aim at developing 
systematic procedures for generating uniformly distributed samples in a given semialge¬ 
braic set. This is an important problem since many system specifications lead to sets with 
a (complex) closed-form description, and being able to draw samples from these type of 
sets provides the means for the design of systems with a complex set of specifications. 
In particular, the algorithm presented in this section can be used to generate uniform 
samples in the solution set of LMIs. 

As before, we assume that the set of interest is a compact basic semialgebraic set defined 
as in ([^, and that there exists a bounding hyper-rectangle B = [a,b] of the form (|^. 
Then, the problem we discuss in this section is the following. 
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Figure 7: Including the same 10 space points (blue) in PSS of degree 4 (left) and 10 
(right). 


Problem 3 (Uniform Sample Generation over /C) Given a semialgebraic set K, de¬ 
fined in a of nonzero volume, generate N independent identically distributed (i.i.d.) 
random samples ..., uniformly distributed in )C. 


Let us start by describing the approach proposed to solve this problem. First, we dehne 
the uniform density over the set K, as follows 


Uyc = 


I/C 

vol 1C 


(16) 


where Ix; is the indicator function of the set JC defined in Q. Then, the idea at the 
basis of the proposed method is to use a PSS approximation of the set /C or, equivalently, 
a polynomial over approximation of the indicator function lyc, obtained employing the 
framework introduced in Sections 2 and 3. 


To this end, given a degree d G N, consider the optimization problem ([^ and let p*^ be a 
polynomial that achieves the optimum. If one examines the proof of Theorem one can 
see that this polynomial has the following properties 


i) P*d > I/c on S 

ii) As d ^ oo, > Ix: both in and almost uniformly on B. 

Hence, p*^ can arbitrarily approximate (from above) the indicator function of the set 1C., 
and therefore it represents a so-called “dominating density” of the uniform density Ux; on 
B. More formally, there exists a value /? > 0 such that fip*d{x) > Ux;((r) for all x £ B. 
Hence, the rejection method from a dominating density, discussed for instance in [HI 
Section 14.3.1], can be applied leading to the following random sampling procedure. 

A graphical interpretation of the algorithm is provided in Figure]^ for the case of a simple 
one-dimensional set 

)C = {xeR : (x - 1)2 - 0.5 > 0, a: - 3 < 0} . 
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Algorithm 1 Uniform Sample Generation in Semialgebraic Set fC 
Given d &N, let be a solution of 

min / p(x)dx 

P^Pd Jb 

s.t. p > 1 on /C 
p > 0 on 15. 

1 . Generate a random sample ^ with density proportional to p^ over B. 

2 . If ^ ^ /C go to step 1. 

3. Generate a sample u uniform on [0, 1]. 

4. If MPd(^) < 1 return x = ^, else go to step 1. 


(17) 


First, problem Q is solved (for d = 8 and B = [1.5, 4]), yielding the optimal solution 

pI{x) = 0.069473x®-2.0515a)^+23.434a)®-139.5x^+477.92a)^-961.88a)3+1090.8x2-606.07x+107.28. 

As it can be seen in Figure p^ is “dominating” the indicator function Ifc on B. Then, 
uniform random samples are drawn in the hypograph of p^. This is done by generating 
uniform samples ^ distributed according to a probability density function (pdf) propor¬ 
tional to p^ (step 2), and then selecting its vertical coordinate uniformly in the interval 
[0, (step 3). Finally, if this sample falls below the indicator function (blue dots) it 
is accepted, otherwise it is rejected (red dots) and the process starts again. 



Figure 8: Illustration of the behavior of Algorithm 1 in the one-dimensional case. Blue 
dots are accepted samples, red dots are rejected samples. 

It is intuitive that this algorithm should outperform classical rejection from the bounding 
set B, since more importance is given to the samples inside K, through the function p^. 
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To formally analyze the performance of Algorithm 1 , we define the acceptance rate (see 
e.g. |E|) as the reciprocal of the expected number of samples that have to be drawn 
from p*^ in order to find one “good” sample, that is a sample uniformly distributed in /C. 
Then, the following result, which is the main theoretical result of this section, provides 
the acceptance rate of the proposed algorithm. 


Theorem 4 Algorithm 1 returns a sample uniformly distributed in K. Moreover, the 
acceptance rate of the algorithm is given by 


yo\K 


Id = 


Wa 


where w*^ = Jjgp'^{x)dx is the optimal solution of problem 


Proof: To prove the statement, we first note that polynomial defines a density 


/ = 


Pd 

w* 


(18) 


over B. Moreover, by construction, we have > Ifc on B, and hence 

A y lit 

vol K, ~ vol /C 

/ ^ 

> —- > 7dU/c 


(19) 


vol /C 


w: 


on B. Then, it can be immediately seen that Algorithm 1 is a restatement of the classical 
Von Neumann rejection algorithm, see e.g. im Algorithm 14.2], whose acceptance rate is 
given by the value of such that (19) holds, see for instance dg. □ 


It follows that the efficiency of the random sample generation increases as d increases, 
and becomes optimal as d goes to infinity, as reported in the next corollary. 

Corollary 2 In Algorithm 1, the acceptanee rate tends to one when inereasing the degree 
of the polynomial approximation, i. e. 


lim = 1 . 

d—)-oo 


Therefore, a trade-off exists between the complexity of computing a good approximation 
{d large) on the one hand, and having to wait a long time to get a “good” sample (7 
large), on the other hand. Note, however, that the first step can be computed off-line for 
a given set /C, and then the corresponding polynomial p^ can be used for efficient on-line 
sample generation. Finally, we highlight that, in order to apply Algorithm 1 in an efficient 
way (step 2 ), a computationally efficient scheme for generating random samples according 
to a polynomial density is required. This is discussed next. 
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8.1 Sample generation from a polynomial density 


To generate a random sample according to the multivariate polynomial density / defined 
in (18), one can use the so-called conditional density method described in [13]. This 
is a recursive method in which the individual entries of the multivariate samples are 
generated according to their conditional probability density. We now elaborate on this. 
We should note that the approach developed in this paper only provides the density up 
to a multiplying constant. However, to simplify the exposition to follow, we proceed as if 
the polynomial given is indeed a probability density function. 

Assume that the bounding set is a hyperrectangle B = [o, b] of the form (|^ and that 
have a polynomial density p. We start by computing the marginal density 


we 


Pi Xi 


rh2 rh 

J a2 J dn 


p{Xi,X2, . . .,Xn) dX2 • • - dXr, 


and, for each i = 2,... ,n and given xi,... Xj_i, compute conditional marginal densities 

dbi^i 


Pi \ Xi ^ 


roi-\-i ron 

J ai-\-i J a-n 


p{x\^ . . . , Xi—\^ Xi^ • • • 5 dXiJ^i • • • dXr^^ 


and respective (polynomial) cumulative distributions F* satisfying 


dxi 

The sampling procedure then starts by computing a sample xi according to Fi and, 
iteratively, computing samples Xi given Xi,.. .Xj^i according to the distribution F). The 
exact description of this procedure is described in Algorithm One should note that, 
given the density p, a closed form is available for all marginal and conditional densities. In 
other words, none of the integrations mentioned above needs to be computed numerically. 


8.2 Numerical example: sampling in a nonconvex semialgebraic 
set 

To demonstrate the behavior of Algorithms 1 and 2, we revisit Example |6.1[ and generate 
uniform samples in the semialgebraic set /C defined in ( [I^ . As already shown in Figure 
|3l the indicator function is well approximated from above by the optimal PSS p*^ ^ for 
a = 20. The results of Algorithm 1 are reported in Figure]^ The red points represent the 
points which have been discarded. To this regard, it is important to notice that also some 
point falling inside /C has been rejected. This is fundamental to guarantee uniformity of 
the discarded points. 


9 Concluding Remarks 

In this paper we have introduced the concept of polynomial superlevel sets (PSS) as a 
tool to construct “simple” approximations of complex semialgebraic sets. Algorithms are 
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Algorithm 2 Generation from a polynomial density 

Returns a sample in with density proportional to the polynomial 


UoL n 

p: Xi,...,Xn ^ n 

i=i i=i 


1. Let i = 1 

2. Compute the univariate polynomial 


ng 

F : Xi ^ '^-fijixi,... ,Xi_i)x'^^’"^^ 

where 


IfiJ (^1 ? • • • 5 l) 


a 




+ 1 


Pj 




1 


ttj-y + 1 



3. Generate a random variable w uniform on [F{ai), F(bi)] 

4. Compute the unique root in [a*, 6*] of the polynomial Xj F{xi) — w 

5. Let Xi = G 

6. If i < n let i = z + 1 and go to (2) 

7. Return x 


( 20 ) 


( 21 ) 


( 22 ) 


provided for computing these approximations. Moreover, it is shown how this concept can 
be used to solve two important problems: i) reconstruction/approximation of sets from 
samples and ii) generation of uniform samples in basic semialgebraic sets. Examples of the 
application of these ideas to problems in control engineering are also described. Note that 
the methods provided in this paper can be used to obtain probabilistic approximations 
of difficult sets, in the spirit of what is discussed in [Hj. Also, in [13] the application 
of minimum size PSS to the approximation of the one-step reachable set of a nonlinear 
discrete-time function is presented, with an extension to nonlinear set filtering. Finally, 
we note that similar techniques can also be used to approximate transcendental (i.e. 
non-semi-algebraic) sets arising in systems control, e.g. regions of attraction, maximum 
positively invariant sets, and controllability regions. 
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